Structure of vortices in rotating Bose-Einstein condensates 
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We calculate the structure of individual vortices in rotating Bose-Einstein condensates in a trans- 
verse harmonic trap. Making a Wigner-Seitz approximation for the unit cell of the vortex lattice, 
we derive the Gross-Pitaevskii equation for the condensate wave function in each cell of the lat- 
tice, including effects of varying coarse grained density. We calculate the Abrikosov parameter, the 
fractional core area, and the energy of individual cells. 
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I. INTRODUCTION 



Rapid rotation of an atomic Bose-Einstein condensate 
has made it possible to create vortex lattices in which the 
vortex core size is comparable to the intervortex spacing 
[Hi S @]- This has opened up a new regime inaccessi- 
ble in superfluid He-II. At low rotation rates, the vortex 
core size is small compared with the intervortex spacing, 
and it is given by the Gross-Pitaevskii healing length 
[3,l5|,lS]- However, as the angular velocity Vl is increased, 
the vortex core radii grow and eventually become com- 
parable to the separation of vortices 0, H, S In 
a harmonic trap, when f2 is close to the transverse trap- 
ping frequency uj^_ , the condensate wave function is dom- 
inated by the lowest Landau level (LLL) of the Coriolis 
force [ll|] . The criterion to be in the LLL regime is that 
TiU, ^ gn, where gn is the interaction energy, n is the 
particle density, g = inh^as/m is the two-body interac- 
tion strength, m is the particle mass, and as is the s-wave 
scattering length. If the system at rest can be described 
by the Gross-Pitaevskii (GP) equation for the condensate 
wave function, then the GP approach is valid for any ro- 
tation rate until the number of vortices becomes of 
order 10 % of the number of particles iV [H [H, Hi, [Hi . 

The purpose of this paper is to calculate the internal 
structure of the vortices in a lattice for arbitrary rotation 
velocities and position in the lattice. Such calculations 
are needed to analyze current experiments in which gn 
is comparable to hfl, and one is between the limits of 
slow rotation and the LLL regime. This work general- 
izes an earlier account Q, which assumed that the vor- 
tex structure was the same throughout the lattice. It 
is also com plem entary to the work of Cozzini, Stringari, 
and Tozzo [iq . who investigated numerically the vortex 
structure in the case when the rotation rate equals the 
transverse trapping frequency, and thus the particle_den- 
sity is uniform. Our basic approach, as in Ref. 8], is 
to treat a unit cell of the vortex lattice in the Wigner- 
Seitz approximation, in which one replaces the actual 
unit cell, hexagonal for a triangular lattice, by a circu- 
lar one. In Sec. [TTl we present the basic formalism. We 
find that the structure depends not only on the rotation 
rate and local density, but also on derivatives of the lo- 



cal density (an effect taken into account only globally in 
Ref. [1]). Numerical calculations of the local structure 
of vortices are given in Sec. IIIII In Sec. IV we calcu- 
late the Abrikosov parameter, fractional core area, and 
the energy of a Wigner-Seitz cell, and discuss the spatial 
dependence of the vortex structure. 



II. BASIC FORMALISM 

In this section, we derive the GP equation for a sin- 
gle Wigner-Seitz cell of the vortex lattice. We assume a 
separable trapping potential. 



V{v) = V^{r^) + V,{z), 



(1) 



where is the coordinate perpendicular to the rota- 
tion axis, which we take to be in the z direction. We 
also assume that the confinement in the z direction 
is sufficiently tight that the system is essentially two- 
dimensional, and that the condensate wave function '5 is 
separable. 



(2) 



We normalize ^ to the number of particles by taking 
J dr±\ip\'^ = N and / dz\H\'^ = 1. In the remainder of 
this paper we consider only a two-dimensional system 
and drop the _L subscript. 

The structure of the condensate wave function in 
the transverse direction is determined by the two- 
dimensional GP equation. 



^V^ + V + g2Bm']i'^Iii^. (3) 
2m ' 



Here fi is the chemical potential, excluding the contribu- 
tions from the motion in the z direction. The effective 
two-dimensional interaction strength is 



92D = g 



dz\H\ 



(4) 



For a system uniform in the z direction, g2T) = g/Z, 
where Z is the extent of the cloud in the z direction, while 



2 



if H is the ground-state wave function for an oscillator 
potential with frequency uj^, g2u = g/V^ndz, where = 
[13. 

We motivate our derivation of the basic equation for 
the vortex structure by focusing on the simple example of 
a vortex at the center of the trap. To solve the GP equa- 
tion, subject to the boundary conditions imposed by the 
presence of the other vortices, we make a Wigner-Seitz 
approximation to the central cell of the vortex lattice, 
replacing the hexagonal unit cell by a circle of the same 
area. The radius of the Wigner-Seitz cell is 



(5) 




WScel 



where is the local vortex density. The wave function 
in the central cell is then a p-wave state, 



(6) 



FIG. 1: A Wigner-Seitz cell of the vortex lattice. The position 
of the center of the cell j relative to the center of the trap is 
Hj , and p is the transverse coordinate measured from Rj . 



where (j) is the azimuthal angle measured relative to the 
center of the trap, and uq is the average of the particle 
density over the central cell. We take the average of |/p 
over the cell to be unity, j\ d^r p — ttP . Without loss 
of generality, we may take / to be real; then / satisfies 



2m 



1_9 

r dr 



,5/ 
dr 



f 



V{r)f 

92Dnaf = (7) 



In order that the wave function in the central cell con- 
nect smoothly with that in neighboring cells, we require 
df /dr = at the boundary of the Wigner-Seitz cell, 
r — £. Note that the structure of the vortex depends 
on the local trapping potential. In the Wigner-Seitz ap- 
proximation, the structure of the vortex on the axis of 
the trap depends on the angular velocity of the trap only 
through the length £ that determines the size of the unit 
cell. If the symmetry of the lattice is taken into account, 
there will be other terms that depend on due to the 
flow created by the rotating lattice of vortices outside 
the central cell, but these will be small because the true 
hexagonal unit cell is well approximated by a circle. 

We now turn to the general vortex lattice. As in Refs. 
0, , we assume that the number of vortices is large and 
write the condensate wave function as the product of the 
square root of a slowly varying coarse-grained density n, 
a rapidly varying real function /, which describes the 
local structure of vortices, and a phase factor e**: 



axis, and that the system is in equilibrium. We then 
specialize, in calculating the vortex structure, to a uni- 
form triangular lattice. In separate publications we will 
investigate the equilibrium distortions of the lattice and 
calculate the elastic constants (18| . 

The local fluid velocity is given by ?iV$/m, which we 
write as a sum of the background flow, v^, plus a locally 
circulating flow: 



h n 

— V$ = i^b + — V(/)j , 



(9) 



where is the azimuthal angle measured with re- 
spect to the center of the vortex in cell j. This 
latter term accounts for the vorticity in cell j. We 
can write the background flow in general as = 
(h/m) [V$(i?j ) -|- V0sj(r)], where Hj is the center of cell 
j and ipsj is regular in the cell [T§| . Here we neglect effects 
due to (l)sj. 

We estimate Vb by calculating the circulation around 
a circle of radius Rj in terms of the number of vor- 



tices Nvj within the circle: 
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hN^j/mR'^. Thus 



ilv^ X Rj, where 



(10) 



(8) 



As in the case of the central cell, we normalize by 
setting its average over each unit cell of the lattice to 
unity. 

We assume, in setting up the general expressions for 
the energy and angular momentum of the vortex lattice, 
that the lattice is close to uniform, with possible displace- 
ments depending only on the distance from the central 



If the lattice is uniform with vortex density , then = 
Ti-KTi^/m. In addition V0j = </>j/p, where (ft^ is the unit 
vector in the azimuthal direction about Rj, and p = 
r — Rj is the coordinate in the x-y plane measured from 
the center of the cell (see Fig. [T|). 

We now derive the Gross-Pitaevskii equation in a 
Wigner-Seitz cell for a rotating condensate in a harmonic 
trap. The total energy in the laboratory frame is [Eq. (8) 
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of Ref. fsj] 



Thus 



E = 



+ — ^ — nf + — n f 



YdV|l|(VV^)V^ + iv/^.Vn| , (11) 



where I = J (Pr ■mn{r)f'^r'^ is the moment of inertia. In 
general, n(r) varies in space, and we continue to allow 
for the possibility that depends on position. 

The angular momentum L [Eq. (11) of Ref. Q] is sim- 
ilarly 



L = J (fr nf[r x hV^], 

- ^J/^^'' "/^N^v(i?j + P ■ Rj) + hp ■ Rj/p^] 
+hN. (12) 



In the remainder of this paper we neglect the weak 
spatial dependence of arising from the small devi- 
ations of the vortex lattice from regular triangularity 
0, HSi im, m, , and assume that the lattice is uni- 
form. We also assume that / is axially symmetric within 
each cell. Our objective is to calculate the energy in a 
frame rotating with angular velocity as a functional of 
/, and then to determine / by minimizing the energy. 
In Eq. pTjl . the energy per particle due to the trapping 
potential {Iuj'^/2) and the kinetic energy of bulk rotation 
(the term) are of order mu'^R^ ~ h{ui'^ /^^)N^ and 
m^lR^ ^ hflvN^, respectively, where R is the radius 
of the cloud. The leading contributions to these terms 
are independent of /, and therefore to determine / it 
is necessary to include the leading /-dependent contri- 
butions, which are smaller by a factor 1/A^v The en- 
ergy per particle from interparticle interactions is of or- 
der Ng2D/R^, and since this depends on / it is unnec- 
essary to include corrections of relative order l/N^. To 
calculate the energy in a uniform lattice, we expand the 
coarse-grained density about the center of the Wigner- 
Seitz cell, n(r) ~ n(Rj) -t- p ■ Vn(Rj) + \{p ■ VfriiHj). 



2m 



dp 



(13) 

Here we have retained in the kinetic and trap energies 
terms of order hfl and hoj'^ /flv per particle, and thus the 
last two terms in Eg. ([TT]). which are smaller by a factor 
1/Ny, are neglected |24l|. To obtain Eq. we used the 
fact that 

nfp ■ r 

(^n{R,) + ii?^.^n(i?,)) J^d'r f p\ 

(14) 

where the factor 1/2 comes from averaging in two dimen- 
sions [2^. 

The angular momentum L [Eq. (fT2|l ] is 
L = j d^r nf[v X 

~ m^ + fiN n{Rj) I d^r mil^fp^ 

-\T.^i^<R^) j <^^'- Mvp'-?^)/'- (15) 



We write the moment of inertia as 
I = 1 + j d'^r mn{f - l)P 



I + mY,[n{Rj)+Rj^n{R 



7 ^ 



V fp' 
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X / d'r if - l)p' 



(16) 



where I = J d^r 'mn{r)r'^ is the moment of inertia calcu- 
lated for the coarse-grained density. Note that here the 
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second derivative term gives contributions of order HQ to 
the energy per particle. In the Umit of a smaU core, i.e., 
for hQ <C g2D'n, I reduces to /. The angular momentum 
written in terms of / is 



L = 



(17) 



where 



and 



niRj) 



n{Rj) dR 



i?2 



d In Ri 



92 



n{R,). 



(18) 



(19) 



From Eqs. (HS]) and ((171), total energy E' in the 
frame rotating with angular velocity is 

E' =E - VlL 



1 



4rn V ri?, 



1 E^.^"(^.) 



(20) 



where we have collected all terms dependent on / in 



+ -mujpf +— n{Rj)f 



(21) 



The effective trap frequency uj is given by 

3 

Co' =uj'' + a{u? - 0.0.^) + ^(cj2 + f^v - Zr^riv)- (22) 

[The coefficient of the term p' p differs from that in Eq. 
(15) of Ref. [1] because of the assumption made there 
that / is independent of position.] The explicit form of 
E is found by setting = in Eq. ([20| . The term pro- 
portional to /(w^ — n^) is the energy due to the trapping 
and centrifugal potentials, while the term proportional 
to /(ri — ilv)^ is the extra kinetic energy due to the fact 
that when 7^ fiv there is bulk motion in the rotat- 
ing frame. The equilibrium value of fiv is determined by 



minimizing E' with respect to flv keeping n and / fixed, 
dE' /dVl^ = 0. This condition has the form 



0, 



(23) 



where the last term is given by E'^.^^ — E' — 
(1/2) [{uj' - n') + {n- fl^)'] and is of order Nhn. 
From this it follows that (£i — ) / fi oc l/N^, in agree- 
ment with earlier works (20l . [21I . [2^ . Since we assume 
that A^v 3> 1, we shall henceforth put ily = fi. Thus 



(24) 



In the fast rotation regime, in which fi ~ cj, the cor- 
rection to the trap frequency due to the a and /3 terms 
is small and a) ~ w for all cells. Thus the vortex struc- 
ture does not change significantly from cell to cell. How- 
ever, at lower rotation velocities, the position-dependent 
a and /3 corrections and the variation of the vortex struc- 
ture from cell to cell are non-negligible. The frequency uj 
equals Lo in the central cell and decreases with increasing 
Rj. 

As explained in earlier works, for Na^/Z ^ 1 — Q,/lo 
d, [21I, [22, Hij (which is equivalent to the condition » 
1) the global density profile is to a good approximation 
a Thomas- Fermi parabola, n(r) = n(0)(l — rfjR') [27|. 
Here R is the radius of the cloud [1, 0, S [ll, and the 
central density is [13, HI] 



„(o) - ^'^ L _ ^ ^'^ 

Here b is the Abrikosov parameter, 

- (/2)2 /. d?r ■ 



(25) 



(26) 



In deriving the Thomas- Fermi form it has been assumed 
that variations of b over the cloud, which amount to at 
the very most 16%, may be neglected. The radius of the 
cloud is given by 



-1/2 



with 



(27) 



(28) 



In the experiments of Ref. [13], = 27r x 5.3 Hz, and 
fls = 5.6 nm for the triplet state of ^^Rb. The number of 
atoms varies from N ~ 10^ for 57 < 0.95w to TV ~ 10^ for 
ri K, 0.99w, and the corresponding values of kq = K'i^/\/b 
are of order 100 and 10, respectively. We calculate below 
for kq = 10 and 100. 
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For a Thomas- Fermi profile, a = f3 
2 [1 - n{0)/n{r)] = -2r^ / [R^ - r^), and [H] 



= 0,2 _ 



2(i?2 



(29) 



Varying Ej with respect to /, subject to the normaliza- 
tion condition ^. d?r p — tt£'^, we obtain the GP equa- 
tion in a Wigner-Seitz cell: 



2m 



1 d 



Of 



pdpydpj p2^ 



/ 



(30) 

This equation, with boundary conditions /(O) = and 
df{£)/dp — 0, determines the optimized form of /. It 
is a generalization of Eq. ([7]) to vortices away from the 
center of the trap. The equation fails for vortices close 
to the edge of the cloud where the change of the coarse- 
grained density within the cell is large and the expansion 
of n(r) about Rj is inadequate; the expansion is justified 
provided 



n{r) 



dn{r) 



dr 



1 



1 



/]Vvi?l-r2/i?2 



« 1 , 



(31) 



a condition that breaks down for the outermost single 
layer of vortices, where r — > i?. 



III. VORTEX STRUCTURE 

We now investigate the vortex core structure using the 
framework derived in the previous section. It is conve- 
nient to measure lengths in units of the effective oscillator 
length d= .Jh/mCj [cf. Eq. jM])]; the GP equation jSO]) 
becomes 



r] dr] \ dr] 



+ ^+ r]\f + - A/ , (32) 



3 
2 

^' 

-1 

-2 



' 'f(i) 

//'"^■x^ f(0) 

I / >r\-^-^' 4°) (central cell) 

/// \-VV^-.. 




FIG. 2: The first five basis functions ff''^ of the free HamiL 
tonian _ffo, for Q, = lo for the central cell, where Cj = uj. 



where Ai is the normalization constant ensuring 
Jj (Pr] [(/i^"*)^ — 1] =0, and iFi is Kummer's conflu- 
ent hypergeometric function (see, e.g., Ref. [3l|). The 
eigenvalue is determined by the boundary condition 



(36) 



where ryccii = £/d = ^Juj/il. Figure [2] shows an example 
of the first five ff*^ plotted as functions of p/£. The 
ith eigenfunction ff^ has i — 1 nodes. At = w and 
thus u) = w, the eigenvalue A^'^-' of the first basis state 



/} , determined by Eq. ([36]), equals 4, and /} reduces 
to the p-wave function for a harmonic oscillator: ff*^ ~ 

(p/c?) exp (— p^/2(i^). Therefore, in this hmit, /^"^ is the 
LLL component of the wave function in the Wigner-Seitz 
cell, and the f'f \ for i>2, describe the higher Landau 
level components. 
We expand / in Eq. 



as 



(0) 



(37) 



where r] = p/d, A = 2mpd'^ /fi^ , and k = 2mg2T:ind^ / . 
The scaled structure of a vortex depends only on the local 
parameters k and (via the boundary conditions) i/d. For 
a Thomas-Fermi density profile. 



K = Kq 



1/2 



1 - 



i?2 



(33) 



We solve Eq. ([32]) by constructing a basis set of eigen- 



functions {fi^^} of the free {k = 0) Hamiltonian Hq — 
-r]^^dy^{rjd^) + l/r;^ + 77^ (see, e.g., Ref. (s^l), satisfying 



(34) 



with eigenvalue The eigenfunction ff'^ of i/o is 



p(0) _ 



A, iFi(l-AfV4,2;r,^ 



776 



-„V2 



(35) 



and obtain the expansion coefficients oi for given k by 
propagating the solution numerically in imaginary time, 
thus damping the contribution of higher energy states. In 
this analysis we use up to 15 basis functions for kq = 10 
and up to 20 basis functions for ko — 100. In general the 
number of basis functions required increases as the core 
size decreases relative to the cell size. 

Before showing results, we comment on the limits of 
the present analysis. The convergence of the expansion 
P7p becomes poor in the slow rotation regime at larger k, 
where higher energy states have significant amplitudes. 
More importantly, we have implicitly assumed [e.g., in 
Eq. ([5])] that the cloud contains a large number of vor- 
tices, which is not the case for slow rotation. From Eq. 
(I27t one finds 



(up' 



-1/2 



(38) 



6 



0.5 



1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 

K„=100 


1 














"y^x-^^ /c.= 10, 100 




1 




(central cell) 








/ n = 0.99cj 




0.5 


- 'A 


n 0.9oj 






- i^'//7 


n = 0.75cj 






ill 

r, 1 , 1 , 


n=o.5w 

, , 1 , , , 1 < , < 1 < , , 


, 






0.2 



0.4 



0.6 

p/L 



0.8 



FIG. 3: (Color online) Vortex core structure at various f2 for 
the central cell obtained from the Gross-Pitaevskii equation in 
a Wigner-Seitz cell, with = 10 (black lines) and ko = 100 
(blue lines). 



which is large compared with unity for fi/cj ^ 

In the following, we illustrate the procedure by cal- 
culating the vortex structure close to the center of the 
trap, where a and (5 can be neglected and Co = u). Figure 
[3] shows the function / obtained from the GP equation 
dSH) in a Wigner-Seitz cell for 0.5 < f^/w < 0.99. With 
increasing fJ, the vortex core radius increases, becoming 
comparable to the Wigner-Seitz cell radius at > COlj. 

In Fig. [3] we compare / at kq = 10 obtained numeri- 
cally with analytic expressions in the fast and slow rota- 
tion regimes. In this figure 



/lll(p) 



1 



(39) 



is the p-wave function of a two-dimensional harmonic os- 
cillator with d = £, the lowest energy solution of Eq. (p2|) 
without the interaction term in the limit ^ —^ uj. 

For slow rotation, the core structure is well approxi- 
mated by the single vortex form Q, 



fcip) = A, 



where the coherence length ^ is given by 



2mg2B n(0) 



and 



l + 2Cln 



2C 



1 + 2C 



-1/2 



(40) 



(41) 



(42) 



with C = ^l^'i normalizes /c according to 

As Fig. IHa) shows, the numerical solution of / at = 
0.99w is well described by /lll- However, the numerical 
solution starts to deviate from /lll with decreasing 
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FIG. 4: (Color online) Comparison of the numerical solution 
of the local wave function / with the LLL structure /lll and 
with the single vortex form /c, calculated in the central cell. 



and, at O = O.Olj, it is better approximated by /c than 
by /lll- For slow rotation the single vortex form /c is a 
good approximation [Figs. Hfc) and|4jd)]. 

Next we investigate the amplitude Ci of the ith basis 
state jf^ in /, Eq. ([57]) . As Fig. O a plot of q for each 
of the cases in Figs. |3] and [4] vs i, shows, c; decreases 
with i almost exponentially, ~ exp [— c(i — 1)], where c 
is determined by kq and At f2 = 0.9w, where the 
system is not yet in the LLL regime, the amplitudes for 
I > 2, corresponding to higher Landau levels (LL) in the 
rapidly rotating limit, are less than 5% (c2 — 0.044); even 
at = 0.5w, they are less than 10% [ci ~ 0.090), so that 
the population in higher LL is less than 1%. A small 
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level index i 

FIG. 5: The amplitudes d of the states f^'^'' in the expansion 
of /, for the cases of acq ~ 10 shown in Figs. [3] and ID In all 
cases here, ci > 0.995 and Ci < 0.1 for i > 2. 
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TABLE I: Abrikosov parameter b and two measures of the 



admixture of order 1% of the higher LL components is 
sufficient to change the internal vortex structure from the 
LLL to the form for an isolated vortex. 



IV. EQUILIBRIUM PROPERTIES OF 
VORTICES 



We now calculate the Abrikosov parameter b, the frac- 
tional area of a vortex core, and the energy of a cell of the 
vortex lattice, focusing on the central cell. In the final 
part of this section we discuss the spatial dependence of 
the vortex structure. 



A. Abrikosov parameter 



Figure [S] shows the Abrikosov parameter as a function 
of Flll = 2nn/g2Dn{0), while TableU gives specific val- 
ues of b for the O values in Figs. [3] and [4] In the small 
core limit, where T^ll ^ 0, 6 = 1. Extrapolation of 
the curve in Fig. [S] to T^^j t = 0, indicates a linear in- 



crease in b at small F^ll, in agreement with Ref. 
The curve increases monotonically with F^ll towards 
b - I.d'r fLjJ^d\ = (e2-5)/[4(e-2)2] 1.158 

as fl —> uj and thus F^ll ~* This value of b is re- 
markably close to that for the exact LLL wave function 
for a triangular lattice, b = 1.1596 [33 |. 



B. Fractional core area 

Reference Q introduced the measure of the fraction of 
the area of a Wigner-Seitz cell taken up by the vortex 
core: 



(43) 



fractional area of the core, r^/£ and /I , at kq 
calculated for the wave functions / in Figs. [3] and 



10, 



Q/lu 


F"^ 

LLL 
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1.00 
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1.158 


0.225 


0.368 


0.99 


2.807 


1.141 


0.217 


0.316 


0.90 


0.826 


1.115 


0.199 


0.231 


0.75 


0.454 


1.087 


0.177 


0.171 


0.50 


0.231 


1.056 


0.117 


0.103 



The factor makes this quantity sensitive to the be- 
havior of f{p) in the outer parts of the cell; in the slow 
rotation regime, the slow variation of / in the outer part 
of the cell contributes significantly to r1/P . In addition, 
for slow rotation (for Vl^QAuj in the present calculations) 
/ develops a local maximum for p < ^, which leads to a 
rapid decrease of rl/l"^ with decreasing Vl at F^ll ^0-2 
for Ko = 10 and at T^l^ < 0.05 for kq = 100. The mea- 
sure ([^5]) of fractional area is also relatively insensitive to 
the variation of the vortex wave function / close to the 
center of the cell. 

An alternative measure of the core size which more 
accurately characterizes the structure close to the vortex 
is 



fii) 



(44) 



In addition, this quantity is not influenced by the behav- 
ior of / in the outer parts of a lattice cell. 

Figure [7] shows these two measures of the fractional 
area of the core for kq = 10 and 100 vs increasing ro- 
tation rate. For slow rotation, the measure more ac- 
curately describes the core size extracted from the data 
[l3| • In the rapidly rotating regime, where the vortex 
core radius is comparable to the intervortex distance, the 
rms radius of the density deficit Tc is readily accessible 
to experiment. For F^^l ^ 0.3, the calculated values of 
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FIG. 7: (Color online) Fractional area of the core for the 
central cell, measured by r^/i?^ (solid lines) and r'^ /i^ (dashed 
lines) (see text), as functions of F^ll. The lines a and c are 
for Ko = 100 and the lines b and d for kq = 10. The dotted 
lines show the empirical expression for the fractional area, 
1.34Fllli obtained for low f2 ^ and the asymptotic value of 
rl/f for large F'^^ i- 



Tc are consistent with the measurements of Refs. 
(see also Table U). As F^ll ~* ^c/^^ approaches 



(cf. r';^/l^ = e- 



_ ll-4e 
~ 6-2e 

0.368 as F 



0.225 



(45) 



1 

LLL 



oo). 



The theoretical estimates of fractional core area are in 
qualitative agreement with the published experimental 
data d, [l^. However, it is important to observe that in 
the analysis of the experiments, core areas were deter- 
mined from a Gaussian fit to the profiles of the density 
deficit in the core. In addition, in our calculations we 
have not allowed for the three-dimensional nature of the 
clouds and the effects of expansion on vortex structure. 
For these reasons, and because of uncertainties in the 
number of particles in the experiments, it is premature 
to perform a detailed comparison between theory and 
experiment. 



C. Energy of the central cell 

In Fig. [HI we show for the central Wigner-Seitz cell, 
the energy Ej=o, the contributions from the kinetic 
and trap energies, i?o,o = Jj^Qd^r n{{h^ /2m)[{dpf)'^ + 
P / p'^]+mu'^ P /2} , and the interaction energy £'o,int = 
(.92d/2) Jj^Qd^f n'^f^, calculated with the numerical so- 
lutions for / as a function of F^ll. The dependence of 
Eofi on r^^L is small, while E'o.int ~ Tlll- Thus in the 
slow rotation limit (r^^L ^ 1), where the interaction 
energy dominates, Eq also roughly scales as Flll- At 
T£^^ ~ 0.54, i?o,o = ^-o.int and £'o,o becomes dominant 
with increasing F^^l- non-interacting limit with 

To exhibit the effect of higher LL components for 
slower rotation, we also show -Eo.b, tbe difference between 
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FIG. 8: (Color online) The energy of the central Wigner- 
Seitz cell Eo (solid line) in units of nh^ /m. The kinetic plus 
trapping potential energy Eo,o (dashed line) and interaction 
energy i?o,int (dash-dotted line) are also plotted. (In the non- 
interacting case, Eq = 2n.) The "binding" energy Eo,^ is the 
energy gain compared to the energy for the first basis state 
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FIG. 9: (Color online) Spatial dependence of the vortex core 
structure for kq = 100 and Q, = 0.9a;. 



Eq and the energy calculated for the first basis state ff'^ . 
At Flll ==0-1' -^o.b is 7% of Eq, faUing to - 0.7% at 
Flll = 1 (at ft = 0.9w and kq = 100, Flll = 0.0826). 
For kq — 10, Eq, i?o,Oj ^-o.int, and i?o,b do not differ sig- 
nificantly from their values for kq ~ 100 at Flll ^ ^■ 
For smaller Flll' both i?o,b and Ea^h/Eo decrease with 
decreasing kq. 
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D. Spatial variation of vortex structure 

Although we have looked so far only at the central 
cell, two effects give rise to spatial dependence of the 
vortex structure. The first is that the coarse grained 
density n depends on position; also the effective trap- 
ping frequency uj depends on spatial derivatives of 
n(r). Both n and uj decrease with increasing distance 
from the center of the cloud, and this increases the 
core radius in the outer region of the cloud. Figure [9] 
shows the variation of the vortex core structure with 
increasing distance r from the center. At r ~ 0.75R the 
core radius is ~ 30% greater than for the central cell. 



V. OUTLOOK 

In this paper we have considered the equilibrium struc- 
ture of vortices in two-dimensional situations. In order 
to be able to make detailed comparison with experimen- 
tal data, it is necessary to take into account the three- 



dimensional nature of the system and the expansion of 
the cloud prior to the measurement of vortex structure. 

In the numerical calculation in this paper we have as- 
sumed a uniform triangular lattice. The next step is to 
extend the methods developed here to consider distor- 
tions of the vortex lattice, to calculate the elastic prop- 
erties of the lattice, including the elastic constants Ci 
and C2 [13, [IB, 28, 33, 34], for general rotation rates in 
the presence of a spatially varying coarse grained density. 



Acknowledgments 

We thank H. M. Nilsen, L. M. Jensen, and H. Smith for 
stimulating and helpful discussions. This work was sup- 
ported in part by Grants-in-Aid for Scientific Research 
provided by the Ministry of Education, Culture, Sports, 
Science and Technology through Research Grant No. 14- 
7939, by the Nishina Memorial Foundation, and by the 
JSPS Postdoctoral Program for Research Abroad. This 
research was also supported in part by NSF Grants No. 
PHY03-55014 and No. PHY05-00914. 



[1] M. R. Matthews, B. P. Anderson, P. C. Haljan, D. S. 

Hall, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. 

83, 2498 (1999). 
[2] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dal- 

ibard, Phys. Rev. Lett. 84, 806 (2000). 
[3] J. R. Abo-Shaeer, C. Raman, J. M. Vogels, and W. Ket- 

terle. Science 292, 476 (2001). 
[4] E. P. Gross, Nuovo Cimento 20, 454 (1961). 
[5] L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 40, 646 (1961) [Sov. 

Phys. JETP 13, 451 (1961)]. 
[6] A. L. Fetter, in Lectures in Theoretical Physics, edited by 

K. Mahanthappa and W. E. Brittin (New York, 1969), 

Vol. XIB, p. 351. 
[7] U. R. Fischer and G. Baym, Phys. Rev. Lett. 90, 140402 

(2003) . 

[8] G. Baym and C. J. Pethick, Phys. Rev. A 69, 043619 

(2004) . 

[9] I. Coddington, P. C. Haljan, P. Engels, V. Schweikhard, 
S. Tung, and E. A. Cornell, Phys. Rev. A 70, 063607 
(2004). 

[10] V. Schweikhard, I. Coddington, P. Engels, V. P. Mogen- 
dorff, and E. A. Cornell, Phys. Rev. Lett. 92, 040404 
(2004). 

[11] T.-L. Ho, Phys. Rev. Lett. 87, 060403 (2001). 

[12] N. K. Wilkin, J. M. F. Gunn, and R. A. Smith, Phys. 

Rev. Lett. 80, 2265 (1998). 
[13] J. Sinova, C. B. Hanna, and A. H. MacDonald, Phys. 

Rev. Lett. 89, 030403 (2002). 
[14] N. R. Cooper, N. K. Wilkin, and J. M. F. Gunn, Phys. 

Rev. Lett. 87, 120405 (2001). 
[15] G. Baym, Phys. Rev. A 69, 043618 (2004). 
[16] M. Cozzini, S. Stringari, and C. Tozzo, Phys. Rev. A 73, 

023615 (2006). 

[17] Y. Castin and R. Dum, Eur. Phys. J. D 7, 399 (1999). 
[18] G. Baym, S. A. Giflord, C. J. Pethick, and G. Watanabe, 



Phys. Rev. A 75, 013602 (2007). 
[19] To make contact between the present decomposition and 
the long-wavelength average in the elastohydrodynamics 
of Refs. [15, 28, 33], we note that the long- wavelength 
average, denoted by {...), of V"l? is mv{r)/h, where ?;(r) is 
the long wavelength flow velocity in the lab. The velocity 
VR{r) in a frame rotating at angular velocity Q is given 
by Uij(r) = i;(r) — f2 X r. This is related to the function 
e(r) which gives the smoothed displacement of a vortex 
from its position for a uniform lattice with vortex density 
riv — mQ,/iTh by vr + 2fl x e(r) — {h/m)W ipa, where 4>s 
describes vorticity-free superfluid flow. Thus (V$(r)) = 
{m/h)fl X [r — 2e(r)] -I- \/(f)s{r), from which we identify 
the long- wavelength limit of 4'sjijc) with the potential 
(jiaijc). The curl of this equation implies V x {V'I>(r)) = 
2(m/?!,)r2(l — V ■ e); since the mean density of vortices 
is given by nv(r) = — V • [nj)e(r)] , we see that V x 
{V$(r)) = 27rnv(r), as expected. Note that this result 
may also be obtained by working in terms of a function 
giving the smoothed displacement of a vortex from the 
eqmUbrium lattice, which has a non-uniform density of 
vortices. 

[20] D. E. Sheehy and L. Radzihovsky, Phys. Rev. A 70, 

051602(R) (2004); 70, 063620 (2004). 
[21] G. Watanabe, G. Baym, and C. J. Pethick, Phys. Rev. 

Lett. 93, 190401 (2004). 
[22] N. R. Cooper, S. Komineas, and N. Read, Phys. Rev. A 

70, 033604 (2004). 
[23] A. Aftalion, X. Blanc, and J. Dahbard, Phys. Rev. A 71, 

023611 (2005). 

[24] Each of these terms is of order h^ /2mR^, and thus they 
scale in the same way with R as does the interaction 
energy ~ g2DN/R'^. They are thus negligible compared 
with the interaction energy provided Na^jZ ^ 1. For 
1 > Nas/Z > 1 - n/uj, the term (n^/2m){V^ff 



10 



would lead to a Gaussian density profile if the vortex lat- 
tice were forced to be perfect However, as has been 
demonstrated earlier, when distortion of the vortex lat- 
tice is allowed for, these contributions are largely can- 
celed by contributions to the elastic energy of the vortex 
lattice and the density profile has the Thomas-Fermi form 
provided ^ 1 [21]. Within the contex of the present 
paper, the Thomas-Fermi form for Nus/Z ^ 1 follows 
from the condition 5E' /Sn = /j., where /i is the chemical 
potential in the rotating frame. 
[25] By contrast, the calculation of Ref. took into account 
only the global average of this term, which vanishes when 
the density vanishes at large distances; as a consequence 
the present equations differ from those of e.g., the 
signs of the second and third terms of Eq. (|16|l are dif- 
ferent from those of Eq. (13) of Ref. [gj] for the global 
average. 

[26] J. Sinova, C. B. Hanna, and A. H. MacDonald, Phys. 

Rev. Lett. 90, 120401 (2003). 
[27] In the present formalism the Thomas-Fermi form for 

Nas/Z 3> 1 follows from the condition 5E'/Sn{r) = 



const with E' given by Eq. (|20|l . To establish this result 
within the present framework for 1 > Na^/Z ^ 1 — Q/uj, 
it is necessary to allow for distortion of the vortex lattice. 

[28] G. Baym, Phys. Rev. Lett. 91, 110402 (2003). 

[29] Note that the effective oscillator potential within a cell 
becomes inverted when n(r)/n(0) < 5(1 — x)/{'^ ~ 5x)i 
where x ~ f2^/<^^; however, this inversion, which at the 
rotations of interest occurs only in the outer layers of the 
lattice, has little practical effect on the vortex structure. 

[30] T. Busch, B.-G. Englert, K. Rzgiewski, and M. Wilkens, 
Found. Phys. 28, 549 (1998). 

[31] Handbook of Mathematical Functions with Formu- 
las, Graphs, and Mathematical Tables, edited by M. 
Abramowitz and I. A. Stegun (National Bureau of Stan- 
dards, Washington, DC, 1972), Chap. 13, p. 504. 

[32] W. H. Kleiner, L. M. Roth, and S. H. Antler, Phys. Rev. 
133, A1226 (1964). 

[33] G. Baym and E. Chandler, J. Low Temp. Phys. 50, 57 
(1983). 

[34] E. B. Sonin, Phys. Rev. A 72, 021606(R) (2005). 



